gusucode.com > 阵列信号处理书的源码 > MATALB 程序/12. 宽带信号DOA估计的非相干信号子空间（ISM）MATLAB程序/ISM_code.m

    %宽带信号DOA估计算法：非相干信号子空间（ISM）算法 
% Developed by xiaofei zhang (南京航空航天大学 电子工程系 张小飞）
% EMAIL:zhangxiaofei@nuaa.edu.cn


clc 
clear all 
close all 
M=12;                                      %阵元数 
N=200;                                    %快拍数 
ts=0.01;                                  %时域采样间隔 
f0=100;                                   %入射信号中心频率  
f1=80;                                    %入射信号最低频率 
f2=120;                                   %入射信号最高频率 
c=1500;                                   %声速 
lambda=c/f0;                              %波长 
d=lambda/2;                               %阵元间距 
SNR=15;                                   %信噪比 
b=pi/180;   
theat1=30*b;                               %入射信号波束角1 
theat2=0*b;                                %入射信号波束角2 
n=ts:ts:N*ts;                             
theat=[theat1 theat2]'; 
 
%%%%%%%%%%%%%%%% produce signal %%%%%%%%%%%%%%%% 
 
s1=chirp(n,80,1,120);                     %生成线性调频信号1； 
sa=fft(s1,2048);                          %进行FFT变换 
%figure, %specgram(s1,256,1E3,256,250);   %频谱图
s2=chirp(n+0.100,80,1,120);                %生成线性调频信号2 
sb=fft(s2,2048);                           %进行FFT变换 
  
% %%%%%%%%%%%%%%%%%%%%% ISM算法 %%%%%%%%%%%%%%%%%% 
P=1:2; 
a=zeros(M,2); 
sump=zeros(1,181); 
for i=1:N
    f=80+(i-1)*1.0; 
    s=[sa(i) sb(i)]'; 
    for m=1:M 
        a(m,P)=exp(-j*2*pi*f*d/c*sin(theat(P))*(m-1))'; 
    end 
    R=a*(s*s')*a'; 
    [em,zm]=eig(R); 
    [zm1,pos1]=max(zm); 
    for l=1:2 
        [zm2,pos2]=max(zm1); 
        zm1(:,pos2)=[]; 
        em(:,pos2)=[]; 
    end 
    k=1; 
    for ii=-90:1:90 
        arfa=sin(ii*b)*d/c; 
        for iii=1:M 
            tao(1,iii)=(iii-1)*arfa; 
        end 
        A=[exp(-j*2*pi*f*tao)]'; 
        p(k)=A'*em*em'*A; 
        k=k+1; 
    end 
    sump=sump+abs(p); 
end 
pmusic=1/33*sump; 
pm=1./pmusic; 
thetaesti=-90:1:90; 
plot(thetaesti,20*log(abs(pm))); 
xlabel('入射角/度'); 
ylabel('空间谱/dB'); 
grid on